#
# More Donors, More Democracy
#
# Author: Sebastian Ziaja
# Date: 27 October 2018
#
# Replication script for
# main text
#
#

workdir <- ""
# workdir <- ifelse(basename(getwd()) == "aiddem_new", "", "../")
source(paste0(workdir, "scripts//load_packages.R"))
source(paste0(workdir, "scripts//extract.custom.R"))
source(paste0(workdir, "scripts//functions.R"))
source(paste0(workdir, "scripts//formulas.R"))
source(paste0(workdir, "scripts//variable_labels.R"))

# determine variable set for core sample generation
vars_core_set <- all.vars(as.formula(bf[[4]]))

sy <- 1994
ey <- 2013

load(paste0(workdir, "data/aiddem_datasets.RData"))
dX <- d1
dd <- read_rds(paste0(workdir, "data/temp/aiddyad.rds"))
dp <- read_rds(paste0(workdir, "data/temp/aidproject.rds"))

ncountries <- dX %>% filter(smpl == 1) %>% select(cc) %>% unique %>% nrow



# Donor proliferation creates a marketplace for idea support ------------


dp %>%
    filter(pur > 14999, pur < 15200) %>%
    group_by(don, rec, year) %>%
    mutate(count = 1) %>%
    summarise(amt = sum(comnet, na.rm = T), prj = sum(count)) -> d2p

d2p %>% ungroup %>%
    mutate(don = 1) %>% group_by(rec, year) %>%
    summarise(prj = sum(prj), don = sum(don)) -> d3p
d3p %>% mutate(avg = prj / don) -> d3p

par(mar = c(4, 4, 1, 1))

# pdf(width = 4, height = 4.6)
boxplot(d3p$avg ~ d3p$don
        , xlab = "Number of donors giving more \nthan 100,000 USD democracy aid"
        , ylab = "Democracy aid projects per donor"
        # , yaxt = "n"
        , ylim = c(0, 17)
        )
# dev.off()



# Empirical analysis -----------------


## Full regression tables for models from the main text --------------

# OLS:

r.ols <- vector("list", 4)
r.ols[[1]] <- felm(v2x.polyarchy.n ~ l.CMgnh |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.ols[[2]] <- felm(v2x.polyarchy.n ~ l.CMgnh +
                    l.CMgap.log + l.CMenh +
                    l.CMeap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.ols[[3]] <- felm(v2x.polyarchy.n ~
                    l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.ols[[4]] <- felm(v2x.polyarchy.n ~ IA_l + l.CMgnh +
                    l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))

screenreg(r.ols
, custom.coef.map =
    list(
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2],
            "IA_l" = vlies[2],
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
    )
            , ci.force = TRUE
        , digits = 2
)

# First stage:

r1 <- vector("list", 4)

r1[[1]] <- felm(l.CMgnh ~ I(l.ZwvCMgwh94 / 10) |
                      cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))
r1[[2]] <- felm(l.CMgnh ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))
r1[[3]] <- felm(I(l.CMgap.log * 10) ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))
r1[[4]] <- felm(IA_l ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))

screenreg(r1
, custom.coef.map =
    list(
            "I(l.ZwvCMgwh94/10)" = vlzp[2],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
    )
        , ci.force = TRUE
        , digits = 2
)

# Second stage:

r2 <- vector("list", 4)

r2[[1]] <- felm(v2x.polyarchy.n ~ 1 |
                      cnamef + periodf |
                      (l.CMgnh ~ l.ZwvCMgwh94) |
                      # (l.CMgnh ~ l.ZwvDistCM) |
                      cnamef, data = dX %>% filter(smpl == 1))
r2[[2]] <- felm(v2x.polyarchy.n ~
                    l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf |
                      (l.CMgnh ~ l.ZwvCMgwh94) |
                      # (l.CMgnh ~ l.ZwvDistCM) |
                      cnamef, data = dX %>% filter(smpl == 1))
r2[[3]] <- felm(v2x.polyarchy.n ~ #1
                    # l.CMgnh +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf |
                      (I(l.CMgap.log * 1) ~ l.ZwvCMgwh94) |
                      cnamef, data = dX %>% filter(smpl == 1))
r2[[4]] <- felm(v2x.polyarchy.n ~ #1
                    # l.CMgnh +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf |
                      (IA_l ~ l.ZwvCMgwh94) |
                      cnamef, data = dX %>% filter(smpl == 1))

screenreg(r2,
       custom.coef.map = list(
            "`l.CMgnh(fit)`" = vlne[2],
            "`I(l.CMgap.log * 1)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
       )
        , ci.force = TRUE
        , digits = 2
)


## Comparison with model without contstitutive terms --------

r1 <- vector("list", 1)

r1[[1]] <- felm(v2x.polyarchy.n ~ l.CMgnh + l.pop.log.r + l.gdpcap.log.r + l.war25
       | cnamef + periodf | 0 | cnamef, data = dX %>% filter(smpl == 1))
r1[[2]] <- felm(v2x.polyarchy.n ~ l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25
       | cnamef + periodf | 0 | cnamef, data = dX %>% filter(smpl == 1))
r1[[3]] <- felm(v2x.polyarchy.n ~ l.CMgnh + l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25
       | cnamef + periodf | 0 | cnamef, data = dX %>% filter(smpl == 1))
r1[[4]] <- felm(v2x.polyarchy.n ~ l.CMgnh : l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25
       | cnamef + periodf | 0 | cnamef, data = dX %>% filter(smpl == 1))
r1[[5]] <- felm(v2x.polyarchy.n ~ l.CMgnh * l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                         cnamef + periodf | 0 | cnamef, data = dX %>% filter(smpl == 1))

(iafe2 <- getfe(r1[[4]])$effect %>% mean %>% sprintf("%2.1f", .))
(iafe1 <- getfe(r1[[5]])$effect %>% mean %>% sprintf("%2.1f", .))



## Distribution of democracy aid amounts by democracy donor proliferation ------

d1 %>% select(CMgnh, CMgap.log) %>%
    filter(CMgnh > 0) %>%
    as.data.frame -> q

par(mar = c(4, 4, 1, 1))

boxplot(q[, 2] ~ q[, 1]
        , xlab = "Number of donors giving more \n than 100,000 USD democracy aid"
        , ylab = "Democracy aid (million USD)"
        , yaxt = "n"
        )


log10.axis(2, at = c(0, 1, 10, 100, 100))
axis(2, labels = c("10", "100", "1,000"), 
     at = c(log10(10), log10(100), log10(1000)))


## Estimated coefficient of the number of democracy donors on democracy -----
## conditional on the amount of democracy aid

rp <- lm(v2x.polyarchy.n ~ l.CMgnh * l.CMgap.log + l.pop.log.r +
                  l.gdpcap.log.r + l.war25 + cnamef + periodf,
         data = dX %>% filter(smpl == 1))

# cluster se at country level

interplot(rp, var1 = "l.CMgnh", var2 = "l.CMgap.log",
          sims = 5000, plot = FALSE, hist = TRUE, xmin = 0) -> ip
# get ratio of non-clustered to clustered coefs:
summary(r1[[5]])$coefficients[1, 2] / (1 * summary(rp)$coefficients[2, 2]) -> iprate
# correct se by ratio:
ip$ub <- ip$coef + ((ip$ub - ip$coef) * iprate)
ip$lb <- ip$coef - ((ip$coef - ip$lb) * iprate)

names(ip)[1:2] <- c("fake", "coef1")  # bug in interplot(); must rename

# plot
interplot(ip, hist = TRUE, var2_dt = dX$l.CMgap.log[dX$smpl == 1]
          ) +
    scale_x_continuous(labels = c("0", "10", "100", "1,000", "")) +
    annotation_logticks(sides = "b") +
    xlab('Democracy aid (m USD)') +
    ylab('Coefficient for the number\n of democracy donors') +
    theme_bw() +
    geom_hline(yintercept = 0, linetype = 2)


# Robustness checks -----------------

# Inidcators of democratic diversity

subx <- c("v2xsharelargestparty", "v2xnopartyplatforms",
          "v2xrespectarguments", "v2xengagedsociety",
          "v2xcivsocpartenviron")

xf <- vector("list", 1)
for (i in 1:length(subx)) {
    xf[[i]] <- sub("v2x.polyarchy.n", paste0("I(", subx[i], " * 10)"), bf[[3]])
}

runmodel(xf) -> r

screenreg(r
          , ci.force = T
          , custom.coef.map = list("`l.CMgnh(fit)`" = vlsets[8])
          , custom.model.names = c("Largest party share", "Party platforms",
                                   "Counterarguments", "Engaged society",
                                   "CSO environment")
)

